!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! Defines a symbolic multivariate monomial.
!!
!! \n
!!
!! In principle, the variables of the monomial can be identified in
!! two ways:
!! \li naming the variables using a corresponding field,
!! \li relying on the positional order.
!!
!! The first option would allow for one variable to appear multiple
!! times into one monomial, but would also imply complications and
!! inefficiencies, so the second option will be used.
!!
!! \note This module should be regarded as part of the polynomial
!! representation provided by mod_sympoly, and not as a stand-alone
!! module. In fact there are many operations that can be done on
!! monomials that produce general polynomials, such as change of
!! variable and integration on the canonical symplex.
!!
!! The main functionalities provided here are:
!! \li constructors: \c symmon, <tt>assignment(=)</tt>
!! \li arithmetic operators: <tt>+, -, *, **</tt>
!! \li point evaluation: \c ev
!! \li derivation: \c pderj
!! \li screen output: \c show
!<----------------------------------------------------------------------
module mod_symmon

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_symmon_constructor, &
   mod_symmon_destructor,  &
   mod_symmon_initialized, &
   t_symmon,       &
   symmon,         &
   assignment(=),  &
   operator(+),    &
   operator(-),    &
   operator(*),    &
   operator(**),   &
   ev,             &
   pderj,          &
   show

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public and private members
 integer, parameter :: &
   wrong_input    = 1, &
   wrong_previous = 2

 !> symbolic monomial
 !! \note \c nsv can be zero, in which case \c degs is allocated to
 !! zero and deg is zero. This is in fact the simplest way to
 !! represent a constant. Alternative representations have
 !! <tt>degs(:).eq.0</tt>, for any size of \c degs.
 !! \note In most of the functions which reference specific variables,
 !! such as evaluation or derivation, it's safe to reference variable
 !! \f$x_j\f$ for a monomial \c m with <tt>j.gt.m\%nsv</tt>, since an
 !! exponent 0 is assumed.
 !! \bug \c nsv should be a derived type parameter
 !<
 type t_symmon
   ! main fields
   real(wp) :: coeff ! coefficient
   integer :: nsv ! number of variables: size of degs
   integer, allocatable :: degs(:) ! partial degrees
   integer :: deg ! = sum(degs)
   ! the following fields should be private on those compilers which
   ! support it
   integer :: ierr = 0 ! error indicator
 end type t_symmon

! Module variables

 ! public members
 logical, protected ::               &
   mod_symmon_initialized = .false.

 character(len=*), parameter :: &
   this_mod_name = 'mod_symmon'

 interface assignment(=)
   module procedure real2symmon
 end interface

 interface operator(+)
   module procedure add
 end interface

 interface operator(-)
   module procedure sub
 end interface

 interface operator(*)
   module procedure mult, mult_scal
 end interface

 interface operator(**)
   module procedure pow
 end interface

 interface ev
   module procedure ev_scal, ev_1d
 end interface

 interface pderj
   module procedure pderj_mon
 end interface

 interface show
   module procedure monshow
 end interface

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_symmon_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_symmon_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_symmon_initialized = .true.
 end subroutine mod_symmon_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_symmon_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_symmon_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_symmon_initialized = .false.
 end subroutine mod_symmon_destructor

!-----------------------------------------------------------------------
 
 pure function symmon(coeff,degs) result(m)
 ! construct a t_symmon object
  integer, intent(in) :: degs(:)
  real(wp), intent(in) :: coeff
  type(t_symmon) :: m

   m%coeff = coeff
   m%nsv = size(degs)
   allocate(m%degs(m%nsv))
   m%degs = degs
   m%deg = sum(m%degs) ! OK: sum is 0 if the size is 0
 
 end function symmon

!-----------------------------------------------------------------------

 elemental subroutine real2symmon(m,x)
  type(t_symmon), intent(out) :: m
  real(wp), intent(in) :: x

   m%coeff = x
   m%nsv = 0
   allocate(m%degs(0))
   m%deg = 0
   
 end subroutine real2symmon
 
!-----------------------------------------------------------------------
 
 elemental function add(m1,m2) result(m)
  type(t_symmon), intent(in) :: m1, m2
  type(t_symmon) :: m
   
   if( (m1%ierr.ne.0).or.(m2%ierr.ne.0) ) then
     m%ierr = wrong_previous
     return
   endif
   if(m1%nsv.ne.m2%nsv) then
     m%ierr = wrong_input
     return
   endif
   if(any(m1%degs.ne.m2%degs)) then
     m%ierr = wrong_input
     return
   endif
   
   m = symmon(m1%coeff+m2%coeff,m1%degs)
 
 end function add
 
!-----------------------------------------------------------------------
 
 elemental function sub(m1,m2) result(m)
  type(t_symmon), intent(in) :: m1, m2
  type(t_symmon) :: m

  real(wp), parameter :: mone = -1.0_wp
 
   m = m1 + mone*m2
 end function sub
 
!-----------------------------------------------------------------------
 
 elemental function mult(m1,m2) result(m)
  type(t_symmon), intent(in) :: m1, m2
  type(t_symmon) :: m
   
  integer, allocatable :: degs(:)

   if( (m1%ierr.ne.0).or.(m2%ierr.ne.0) ) then
     m%ierr = wrong_previous
     return
   endif

   allocate(degs(max(m1%nsv,m2%nsv)))

   if(m1%nsv.ge.m2%nsv) then
     degs(1:m2%nsv) = m1%degs(1:m2%nsv) + m2%degs(1:m2%nsv)
     degs(m2%nsv+1:m1%nsv) = m1%degs(m2%nsv+1:m1%nsv)
   else
     degs(1:m1%nsv) = m1%degs(1:m1%nsv) + m2%degs(1:m1%nsv)
     degs(m1%nsv+1:m2%nsv) = m2%degs(m1%nsv+1:m2%nsv)
   endif
   m = symmon(m1%coeff*m2%coeff,degs)

   deallocate(degs)

 end function mult
 
!-----------------------------------------------------------------------

 elemental function mult_scal(r,m1) result(m)
  real(wp), intent(in) :: r
  type(t_symmon), intent(in) :: m1
  type(t_symmon) :: m

   if(m1%ierr.ne.0) then
     m%ierr = wrong_previous
     return
   endif

   m = symmon(r*m1%coeff,m1%degs)
 
 end function mult_scal
 
!-----------------------------------------------------------------------

 elemental function pow(m1,n) result(m)
  type(t_symmon), intent(in) :: m1
  integer, intent(in) :: n
  type(t_symmon) :: m

   if(m1%ierr.ne.0) then
     m%ierr = wrong_previous
     return
   endif

   m = symmon(m1%coeff**n,n*m1%degs)
 
 end function pow
 
!-----------------------------------------------------------------------

 !> Monomial evaluation. <tt>size(x)</tt> can be bigger than
 !! <tt>m\%nsv</tt>, in which case the missing exponents are treated
 !! as 0.
 !<
 pure function ev_scal(m,x) result(z)
  type(t_symmon), intent(in) :: m
  real(wp), intent(in) :: x(:)
  real(wp) :: z

   z = m%coeff * product(x(1:m%nsv)**(m%degs)) ! OK: size 0 => product 1
 
 end function ev_scal

!-----------------------------------------------------------------------

 !> See \c ev_scal. <tt>size(x,1)</tt> can be bigger than
 !! <tt>m\%nsv</tt>, in which case the missing exponents are treated
 !! as 0.
 !<
 pure function ev_1d(m,x) result(z)
 ! first index: different variables
 ! second index: different points
  type(t_symmon), intent(in) :: m
  real(wp), intent(in) :: x(:,:)
  real(wp) :: z(size(x,2))

  integer :: i
 
   do i=1,size(z)
     z(i) = ev(m,x(:,i))
   enddo
 
 end function ev_1d
 
!-----------------------------------------------------------------------
 
 !> Compute \f$\partial_{x_j}m\f$ for a monomial \f$m\f$. \c j can be
 !! larger than <tt>m\%nsv</tt>, in which case the result is zero.
 !<
 elemental function pderj_mon(m,j) result(djm)
 ! compute the j-th partial derivative
  type(t_symmon) :: djm
  type(t_symmon), intent(in) :: m
  integer, intent(in) :: j

  ! being an elemental function, we can't use degs(m%nsv)
  integer, allocatable :: degs(:)
 
   if(m%ierr.ne.0) then
     djm%ierr = wrong_previous
     return
   endif

   if(j.gt.m%nsv) then
     djm = 0.0_wp
   else
     if(m%degs(j).eq.0) then
       djm = 0.0_wp
     else
       allocate(degs(m%nsv))
       degs = m%degs
       degs(j) = m%degs(j)-1
       djm = symmon(m%coeff*real(m%degs(j),wp),degs)
       deallocate(degs)
     endif
   endif
 
 end function pderj_mon
 
!-----------------------------------------------------------------------

 subroutine monshow(m)
  type(t_symmon), intent(in) :: m
 
  integer :: i, is, ie
  integer, parameter :: clen = 2+1+1
  character(len=9+(clen+2)*m%nsv) :: monchar

   select case(m%ierr)
    case(0)

     monchar = ''
     write(monchar(1:9),'(E9.2)') m%coeff
     ie = 9
     do i=1,m%nsv
       if(m%degs(i).gt.0) then ! don't print x^0 terms
         is = ie+1
         if(m%degs(i).lt.10) then
           ie = is-1 + clen + 1
           write(monchar(is:ie),'(A,I1,A,I1)') ' x',i,'^',m%degs(i)
         else
           ie = is-1 + clen + 2
           write(monchar(is:ie),'(A,I1,A,I2)') ' x',i,'^',m%degs(i)
         endif
       endif
     enddo

     write(*,'(A)') trim(monchar)


    case(wrong_input)
     write(*,*) 'wrong_input'
    case(wrong_previous)
     write(*,*) 'wrong_previous'
   end select

 end subroutine monshow
 
!-----------------------------------------------------------------------
 
end module mod_symmon

